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ABSTRACT 

In this paper we present a framework which provides an analytical (i.e., infinitely dif- 
ferentiable) transformation between spatial coordinates and orbital elements for the 
solution of the gravitational two-body problem. The formalism omits all singular vari- 
ables which otherwise would yield discontinuities. This method is based on two simple 
real functions for which the derivative rules are only required to be known, all other 
applications - e.g., calculating the orbital velocities, obtaining the partial derivatives of 
radial velocity curves with respect to the orbital elements - are thereafter straightfor- 
ward. As it is shown, the presented formalism can be applied to find optimal instants 
for radial velocity measurements in transiting exoplanetary systems to constrain the 
orbital eccentricity as well as to detect secular variations in the eccentricity or in the 
longitude of periastron. 

Key words: Methods: Analytical - Celestial Mechanics - Ephemerides - Stars: 
binaries - Techniques: radial velocities 



1 INTRODUCTION 

In recent years, precise radial velocity (RV) observation of 
stars and careful analysis of RV data have become rele- 
vant in astrophysics since the vast majority of the known 
extrasolar planets have been discovered using this method 
jMavor fc Queloj Il995l ) or have been confirmed this way 
when the planet itself wa s first detected as a transiting ob- 
ject around its host star (|Konacki et aLll2003l ). In the cases 
where there are no observable transits, the characterization 
of extrasolar planets relies on the radial velocity observations 
alon43- Moreover, observations of photometric transits in ad- 
dition to radial velocity measurements constrain the mass 
of the planet (instead of yielding only a lower limit), and 
provide more precise information on the epoch and period; 
see e .g. the case of the low-mass transiting planet HAT-P- 
11b (|Bakos et al.ll2009l ). where the uncertainty in the epoch 
would be ~ 500 times larger if the analysis had relied only 
on radial velocity measurements. Thus, incorporating con- 
straints given by transit timings reduce the uncertainties in 
the RV amplitude and the orbital parameters (e.g., semima- 
jor axis, eccentricity). 

The aim of this paper is to present a set of analytic re- 
lations (based on a few smooth functions defined in a closed 



* E-mail: apal@szofi.net 

^ Of course, with the exception of planets discovered by mi- 
crolensing or direct imaging. The characterization scheme of plan- 
ets around pulsars are highly similar to the analysis of RV data. 



form) which provides a straightforward solution of Kepler's 
problem, and consequently, time series of RV data and RV 
model functions. Due to the analytic property, the partial 
derivatives can also be obtained directly and therefore can 
be utilized in various fitting and data analysis methods, in- 
cluding, for instance, the Fisher analysis of covariances, un- 
certainties and correlations. The functions presented here 
are nearly as simple to manage as trigonometric functions. 

As an application, we give a detailed description about 
scheduling radial velocity measurements in the case of tran- 
siting extrasolar planets in order to derive accurate orbital 
eccentricity and/or detect the variations in the eccentricity. 
Discussions related to this problem in the case of extrasolar 
planets with no additional constraints on their orbital mo- 
tions other than RV d ata can be foun d in lLoredo fc Chernofil 
( 20031) . iFordl (|2008l ) or lBaluevI (|2008l ). Photometric measure- 
ments for transits constrain the orbital period and epoch 
much more precisely than pure RV observations. Thus, for 
a given fixed eccentricity and argument of pericenter, the 
optimal time instances for RV observations depend only on 
the orbital phase. In the case of transiting planets the long- 
term variations in the orbital eccentricity e and longitude of 
pericenter are quite relevant. The presented strategy fo- 
cuses on the precise measurement of the Lagrangian orbital 
elements k = e cos and h = e sin w (including the cases 
when the cicrular property is intended to be confirmed at 
high significance). 

In Section (2] the basics of the mathematical formal- 
ism are presented, including the rules for calculating partial 
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derivatives. In Section O the solution of the spatial problem 
is shown, supplemented with the inverse problem, still using 
infinitely differentiable functions. The spatial case discusses 
how the inclination and the argument of node should be 
incorporated in the formalism without loosing the analytic 
property while the inverse problem describes how the orbital 
elements are derived from the coordinates and velocities. In 
Section |4l we demonstrate how this formalism can be used 
for transiting extrasolar planets to efficiently schedule the 
phase of RV observations in order to minimize uncertainties 
in the Lagrangian orbital elements. The results are summa- 
rized in the last section. 



2 MATHEMATICAL FORMALISM 

The solution for the time evolution of Kepler's problem can 
be derived in th e standard way as given in various text- 
books (see, e.g., iMurrav fc Dermottlll999l ). The restricted 
two body problem itself is an integrable ordinary differen- 
tial equation. In the planar case, three independent integrals 
of motion exist and one variable has uniform monotonicity. 
The integrals are related to the well known orbital elements, 
which are used to characterize the orbit. These are the semi- 
maior axis a, the eccentricity e and the longitude of pericen- 
teo m. The fourth quantity is the mean anomaly M — nt, 
where n = ^fjifc^ = 2-k/ P, the mean motion, which is zero 
at pericenter passag^f] and t is the elapsed time since the 
pericenter passage. The solution to Kepler's problem can be 
given in terms of the mean anomaly M as defined as 



E - esmE = M, 



(1) 



where E is the eccentric anomaly. The planar coordinates 
are 

^ = ifo cos CC7 — r;o sin-cij, (2) 
ri = S^osmvj + rjocosvj, (3) 

where 

Co = a(cos£;-e), (4) 
7)0 = a\/l — sin is ; (5) 

see also iMurrav fc DermottI (|l999l ). Sect. 2.4 for the deriva- 
tion of these equations. Since for circular orbits the longi- 
tude of pericenter and pericenter passage cannot be defined, 
and for nearly circular orbits, these can only be badly con- 
strained; in these cases it is useful to define a new variable, 
the mean longitude as, \ = M -\- uj to use instead of M. 
Since zu is an integral of the motion, A = M = n. Therefore 
for circular orbits zu = Q and equations ©-((Sj should be 
replaced by 



Co = a cos A, 
?7o = a sin A. 



(6) 
(7) 



^ In two dimensions, the argument of pericenter is always equal 
to the longitude of pericenter, i.e. ro = 

^ Throughout this paper the mass parameter of Kepler's problem 
is denoted by /i = Q{mi +m,2), where mi and m2 are the masses 
of the two orbiting bodies and Q is the Newtonian gravitational 
constant. The orbital period is denoted by P. 



To obtain an analytical solution to the problem, i.e. which 
is infinitely differentiable with respect to all of the orbital 
elements and the mean longitude, first let us define the La- 
grangian orbital elements k = ecosro and h — esinzu. Sub- 
stituting equations ©-([S]) into equations (HJ-Q gives 



where c 

1 - x/r= 



c \ e sin E I +h 
s I ^ 2-£ \-k 



(8) 



cos(A -I- esini?), s = sin(A -I- esini?) and £ = 
2, the oblateness of the orbit. The derivation of 
the above equation is straightforward, one should only keep 
in mind that E + zu = X + e sin E. In the first part of this 
section we prove that the quantities 



p{X,k,h) = 



and 



q{X,k,h) 



if fc = and /i = 

e sin E otherwise 



if fc = and /i = 

e cos E otherwise 



(9) 



(10) 



are analytic - infinitely differentiable - functions of A, A; and 
h for all real values of A and for all k^ + = < 1. 
In the following parts, we utilize the partial derivatives of 
these analytic functions to obtain the orbital velocities and 
we derive other relations. In this section we only deal with 
planar orbits, the three dimensional case is discussed in the 
next section. 



2.1 Partial derivatives and the analytic property 

A real function is analytic when all of its partial derivatives 
exist, the partial derivatives are continuous functions and 
only depend on other analytic functions. It is proven in Ap- 
pendix |A] that the partial derivatives of g = q{X,k,h) and 
p = p(A, k, h) axe the following for (fc, h) 7^ (0, 0): 



dq 

dX 
dq 
dk 
dq 
dh 

and 

dp 

dX 

dp 

dk 

dp 

dh 



+s 

1-g 



cos(A + p) ~ k 
sin(A + p) — h 



+ sin(A + p) 
— cos(A -I- p) 



1 



(11) 
(12) 

(13) 

(14) 
(15) 
(16) 



Since for all k^ + < 1, q < 1 and therefore 1 — g > 0, 
all of the above functions are continuous on their domains. 
Since the sin(-) and cos(-) functions are analytic, therefore 
one can conclude that the functions q(-, ■, ■) and p(-, ■, ■) are 
also analytic. 

Substituting the definition of p = p(A, k, h) into equa- 
tion ((8]), one can write 



cos(A + P) ] , P ( +h 
sin(A +p) ) 2-£\-k 



(17) 
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while the radial distance of the orbiting particle from the 
center is \/ + rj^ = r = a(l — g). For small eccentricities in 
equation (|17|) third term (fc, h) is negligible compared to the 
first term (cos, sin) while the second term {h, —k)p/{2 — £) is 
negligible compared to the third term. Therefore for e ^ 1, 
p is proportional to the difference between the true anomaly 
and the mean longitude and q is proportional to the distance 
offset relative to a circular orbit; both caused by the non- 
zero orbital eccentricity. 

Since equation (I17p is a combination of purely analytic 
functions, the solution of Kepler's problem is analytic with 
respect to the orbital elements a, {k,h), and to the mean 
longitude A in the domain a > and + < 1. We note 
here that this formalism omits the parabolic or hyperbolic 
solu tions. The formalism ba sed on the Stumpff functions 
(see IStiefel fc Scheifeld Il97ll l provides a continuous set of 
formulae for the elliptic, parabolic, and hyperbolic orbits 
but this parametrization is still singular in the e ^ limit. 

2.2 Orbital velocities 

Assuming a non-perturbed orbit, i.e. when (fc, h) = 0, and 
d = and when the mean motion n = A is constant, the 
orbital velocities can be directly obtained by calculating the 
partial derivative of equation ((17} with respect to A and 
applying the chain rule since 



d_ U 



dX_d_ 



(18) 



Substituting the partial derivative equation p4p into the 
expansion of d^/d\ and drj/dX one gets 



1-q 



— sin(A -l-p) 
+ cos(A +p) I ' 2- e 



+ 



+h 

-k 



(19) 



Note that equation (|19p is also a combination of purely an- 
alytic functions, the components of the orbital velocity are 
analytic with respect to the orbital elements a, {k,h), and 
to the mean longitude A. 

It is also evident that the time derivative of equa- 
tion ([191) is 



(l-g)« 



cos(A + p) 
sin(A + p) 



+ 



(20) 



+h 

-k 



Obviously, equation (|20|l can be written as 



(l-g)3 U 



(21) 



which is equivalent to the equations of motion since p — 
n^a? and vs+^ — ^ — a(l 



2.3 Other properties 



In this subsection we summarize some other properties of the 
functions p = p(A, k, h) and q = q(A, k, h) which can also be 
helpful in some derivations or during numerical evaluation. 

One of the most important properties is the rotational 
invariance. This is a direct consequence of the relation M — 



X — zij, i.e. q and p would not change if the mean longitude 
is increased by an arbitrary angle of Q and simultaneously 
the vector (fc, h) is rotated with the same angle. Therefore, 

q — q{X — ^l, k mail + hsini}, —ksinQ + h cos Q), (22) 

p — p{X — fl, k cosil + hsinfl, —ksinQ + hcosQ). (23) 

This property also results that 

q — q(0, fc cos A + /i sin A, — fc sin A + /i cos A), (24) 

p — p(0, fc cos A -|- /i sin A, — fc sin A -f /i cos A). (25) 

The terms c = cos(A -I- p) and s = sin(A -I- p) appear 
frequently in the expressions for both the coordinates and 
the velocities. These are related to the {q,p) functions and 
the orbital elements (fc, h) as 

k — qc + ps, (26) 
h — qs ~ pc (27) 
or similarly 

q = kc + hs, (28) 
p = ks- he. (29) 

The partial derivatives of c and s with respect to the mean 
longitude and the Lagrangian orbital elements k and h are 



d_lc 

dX ys 

and 

d{c, s) 
d{k,h) 



1 



1 — q \ — sc 



(30) 



(31) 



3 THREE DIMENSIONAL CASE AND THE 
INVERSE PROBLEM 

For an eccentric and inclined orbit, the spatial coordinates 
of an orbitin g body can be obt ained in the similar manner as 
it is done in iMurrav fc Dermot t (1999), Sect. 2.8. Without 
going into the details, here we present the results of the basic 
calculations. The spatial coordinates r = (a;, y, z) are 



r = P3P2Piro, 



(32) 



where ro = (Coif^OjO) (see also equations |3] and O , and Pi, 
P2 and P3 are the rotational matrices with respect to the 
argument of pericenter, uj, the inclination, i, and the argu- 
ment of the ascending node, n. By substituting the equa- 
tions for {(^o,rjo) into equation (|32p and using the rotational 
invariance of q(-,-,-) and p(-,-,-), i.e. equations (|22p - (|23p . 
one can derive the spatial coordinates of the orbiting body 
in a similar manner as it was done in the planar case. 
The result can be written in a compact form using some 
other auxiliary quantities. First, define the Lagrangian or- 
bital elements ix = 2 sin(i/2) cos SI, iy — 2 sin(i/2) sin fi, 
iz ~ \/4:~ ix ~ iy, and the quantity W — rjix ~ ^iy, where 
{^,ri) is defined in equation ((8)1. The spatial coordinates are 
then 



(33) 
(34) 
(35) 
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Since is a linear combination of ^ and 77, the orbital ve- 
locities are linear combinations of ^ and r), with the same 
coefficients, i.e. 


i = t+]^iyW, 


(36) 


y = n~\ixW, 


(37) 


i = \i^W, 


(38) 


where W = rji^ ^iy 




3.1 The inverse problem 




To compute the orbital elements (a,\,k,h,ix,iy) from the 
spatial coordinates {x,y,z) and velocities {x,y,z), first de- 
fine 


2 2,2,2 
r = X + y + z , 

2 . 2 , . 2 , .2 

V = X + y + z , 

Cx = - zy, 
= xy- yx, 

^2 1 ^2 1 ^2 

c = XX -\- yy ~\- zz. 


(39) 
(40) 
(41) 
(42) 
(43) 
(44) 
(45) 



The Lagrangian orbital elements for the inclination and the 
argument of the ascending node is then 

^ (46) 



V2 Cx 



(47) 



For the eccentricity and the longitude of pericenter one gets 



M \ -I- 



1 X 



(48) 

I r \ II - — : r- _ 1 

The semimajor axis sa tisfies the well known relation 
l|Murrav fc Dermott|[l999l ') 

2\ -1 



2 _ 7/ 

/i(l — fc2 — h"^) \r 11 
The mean longitude is then 



(49) 



A = arg ry 



rzcy 



rzc 



he 



C + c, 2-£j C 



(50) 



where 1=1 — \/l — k"^ — h"^. It can be shown that in the 
planar case, i.e. when z = and i = 0, equation (148 1) and 
equation (|50[) do reduce to 



(51) 



hj ryy 
and 



, . he . kc 
X = arg I ry + ; , —rx 



2-t 



2 - . 



(52) 



respectively. Here, arg(-,-) is defined as arg(a;, y) 
arctan(i//a;) if s ^ and tt -|- arctan(i//a;) otherwise. 




0.4 0.6 
Phase 

Figure 1. Radial velocity variations for a transiting planet or- 
biting its star on a circular orbit. The transit occurs at zero (or 
unity) phase. The four dots represent the phases ip = 0.1292, 
0.4138, 0.5862 and 0.8708 when RV measurements should be ob- 
tained to achieve the largest significance of the orbital circularity. 



4 APPLICATIONS 

The utilization of the presented formalism can cover various 
aspects of RV curve analysis. Due to the conventions and the 
definitions of the orbital elements used in the astrophysics of 
stellar binaries, the evaluation of equation (|19|l simply yields 
the actual radial velocity as the rj component Q. The partial 
derivatives in equations (lll|l - (|16|l and in equations (|30|l - (|31| ) 
can easily be used to calculate the parametric derivatives of 
the RV curves and therefore to support various fitting a l- 
gorithms (e.g. Levenberg-Marquardt, see lPress et al.|[l992l ). 
Moreover, the secular variations in the radial velocities due 
to the secular changes in the orbital elements can also be 
estimated in an analytic way. 

In the follow-up observations of planets discovered by 
transits in photometric data series, the detection of system- 
atic variations in the RV signal is one of the most relevant 
steps, either to rule out transits of late-type dwarf stars, 
and/or blends, or to characterize the mass of the planet 
and the orbital parameters. Since transit timing constrains 
the epoch and orbital period much more precisely than ra- 
dial velocity alone, these two can be assumed to be fixed 
in the analysis of the RV data. However, this constraint 
also includes an ad ditional feature. Using equation (1) in 
iPal fc Kocsii (|2008l ). the mean longitude at transits can be 
calculated as 

k+-^.l + h-J^]-'-^, (53) 



Atr = art 



2-e' 



2-e 



h 



therefore the mean longitude at the orbital phase (p becomes 
A = Atr + 2TTip. Consequently, the partial derivatives of the 
rj RV component v = rj{\ti + 2-nifi, k, h) with respect to the 
orbital elements k and h are 



dv 
dk 
dv 
dh 



drj drj dXti 

dk ~^ dxlm 

drj drj d\ti 
dh ~^ 'd\~dh 



(54) 
(55) 



* In this case an = K/^/l — e'^, where K is the semi-amplitude 
of the RV curve 
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A radial velocity curve of a star, caused by the pertur- 
bation of a single companion can be parametrized by six 
quantities: the semi-amplitude of RV variations, K, the zero 
point, G, the Lagrangian orbital elements, {k, h), the epoch, 
To (or equivalently the phase at an arbitrary fixed time in- 
stant) and the period P. In the cases of transiting planets, 
the later two are known since the transit photometric obser- 
vations constrain both with exceeding precision (relative to 
the precision attainable purely by the RV data). Therefore, 
one has to fit only four quantities, i.e. a — {K, G, k, h). The 
vast majority of the known transiting planets orbits their 
host stars on a tight orbit and these tight orbits are expected 
to be circular (however, there are few known exceptioni|f|) ■ 

Now, we show how the optimal phases of observations 
can be determined to confirm the orbital circularity, i.e. re- 
sult the smallest uncertainty in the orbital elements k and 
h. To obtain the uncertainty of a fitted parameter and/or 
the correlation between t he paramete r s, the Fisher matrix 
method can be utilized (|Finnl 1 19921 : iBaluevI [200^ 1. This 
method gives the covariance matrix as 



{Sa,nSan) ■- 
where 

r - V 



dmf{aL,ti)d„f{a,ti) 



(56) 
(57) 



Here /(a, t) is the model function which depends on its ad- 
justed parameters a = (ai, . . . , ajv), and t represents the 
independent variable(s) (in the case of time series, there is 
one independent variable, the time itself). Since in our case 
there are four unknowns, one has to have at least four data 
points, in order to completely determine the parameters. 
Nevertheless, equation (|57p can be formally evaluated pro- 
viding both the uncertainties and correlations. To determine 
the best four phases (^i, ip2, fa and (^4) which minimizes 
the uncertainty in the eccentricity, one has to minimize the 
volume of the covariance ellipsoid of the (fc, h) p arameters 
(a.k.a. generalized _D-optimality, see lBaluevll2008h . It can be 
shown that this volume is 



U 



dot 



(Sk^) 
(ShSk) 



{SkSh) 



1/2 



(58) 



{(5fc2) - {SkShf 



where the covariances can be calculated using equation (I56|) . 
i.e. (<5fc2) = (r-i)33, (SkSh) = (r-i)34, and {Sh"") = 
(r~^)44 if the parameters are a = (ai, a2, 03, 04) = 
{K,G,k, h). 

We minimized U as the function of the phases ipi , . . . , 
ifii usin g the Mark ov chain Monte-Carlo (MCMC) method 
(see e.g. lFord|[2003 ). Multiple chains were initiated from four 
random phases (namely, all of them are chosen uniformly 
from the interval [0, 1]). We found that all of the chains 
converge to a single set of phases that represent the min- 
imum of the function U(ipi,ip2,'Pz,f4)- Therefore, the set 
of optimal phases is unique: we found the optimal phases 
are ip = 0.1292, 0.4138, 0.5862 and 0.8708. In Fig. [f|these 
phases are marked on a hypothetical radial velocity curve. 
We note that this volume of the covariance ellipsoid of the 



see e.g. http://www.exoplanet.eu for up to date information 



Table 1. Optimal phases of radial velocity measurements in or- 
der to obtain the orbital eccentricity as precise as possible. The 
eccentricity and the alignment of the orbit are quantified by the 
Lagrangian orbital elements k and h. 



k 


h 




•fi2 




ip4 


—0.4 


—0.4 


0.1305 


0.2064 


0.2519 


0.6943 


—0.4 


—0.2 


0.1060 


0.2048 


0.2847 


0.7985 


—0.4 


0.0 


0.0787 


0.1879 


0.3125 


n 8695 


—0.4 


0.2 




0.1584 


3398 


0.9197 


—0.4 


0.4 


0.0316 


0.1180 


0.3701 


n 9555 


—0.2 


—0.4 


0.1964 


0.3307 


0.3910 


0.7027 


—0.2 


—0.2 


0.1497 


0.3180 


0.4207 


0.7943 


—0.2 


0.0 


0.1076 


0.2927 


0.4522 


0.8616 


—0.2 


0.2 


0.0722 


0.2551 


0.4900 


0.9113 


—0.2 


0.4 


0.0437 


0.2040 


5399 

KJ . Kj'Jiy ij 


0.9481 


0.0 


—0.4 


0.2557 


0.4672 


n 5328 


0.7443 


0.0 


—0.2 


0.1854 


0.4445 


n 5555 


0.8146 


0.0 


0.0 


0.1292 


0.4138 


0.5862 


0.8708 


0.0 


0.2 


0.0850 


0.3728 


0.6272 


0.9150 


0.0 


0.4 


0.0511 


0.3169 


0.6831 


0.9489 


0.2 


-0.4 


0.2973 


0.6090 


0.6693 


0.8036 


0.2 


-0.2 


0.2057 


0.5793 


0.6820 


0.8503 


0.2 


0.0 


0.1384 


0.5478 


0.7073 


0.8924 


0.2 


0.2 


0.0886 


0.5100 


0.7449 


0.9278 


0.2 


0.4 


0.0519 


0.4601 


0.7960 


0.9563 


0.4 


-0.4 


0.3057 


0.7481 


0.7936 


0.8695 


0.4 


-0.2 


0.2016 


0.7153 


0.7952 


0.8939 


0.4 


0.0 


0.1305 


0.6875 


0.8121 


0.9212 


0.4 


0.2 


0.0803 


0.6602 


0.8416 


0.9467 


0.4 


0.4 


0.0445 


0.6299 


0.8820 


0.9684 



(fc, h) parameters is approximately 12 times smaller on the 
average if the phases were chosen randomly, and 2.5 times 
smaller if the phases were chosen to be closer with a factor 
of 2 to the phases 0.25 and 0.75 (i.e. <fi' = 0.1896, 0.3319 
and 0.6681, 0.8104, respectively). 

Of course, the same kind of calculation of the phases 
that yields the smallest combined uncertainty in the (fc, h) 
parameters can be performed for arbitrary orbital eccentric- 
ity. For some certain values of {k,h), these optimal phases 
are shown in Table [1] Thus, if some initial values for the or- 
bital elements k and h are known, further observations can 
be planned accordingly. 

Similarly, the optimal phases can be derived for 
arbitrary number of observations (A^obs). However, it 
turns out that for Nobs ^ 5, the phase volume U have 
more than one local minima. In order to find the global 
minimum, we initiated several hundreds or thousands 
of individual initial conditions and applied both the 
previously discussed MCMC method and the downhill 
simplex algorithm (Press et al.l il99^ ) in order to find 
the global minimum of U. For Nohs ~ 5, there are two 
local minima and the global minimum is at the phases 
<fi = (0.1318,0.3978,0.5,0.6022,0.8682). For N^bs > 5, 
one or more of the phases will be degenerated. In other 
words, one should take RV measurements at the same 
phases in order to minimize the uncertainty in the or- 
bital eccentricity. For instance, the global minimum is 
at ip = (0.1376,0.4204,0.4204,0.5796,0.5796,0.8624) 
for Aobs ~ 6 or (y3 = 
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(0.1405, 0.4315, 0.4315, 0.5965, 0.5965, 0.8746, 0.8746) for 
A^'obs = 7. It can be shown that for larger A'^obs which are 
multiple of 4, the optimal phases will be at the exactly same 
location as in the case of A'obs = 4 (see earlier or Table [l|, 
and at each phase one should acquire the same number of 
measurements, namely Afobs/4. 

Due to the limitations in the telescope time and in 
the observation conditions (including day/night variations 
or the visibility of the target object), RV observations can- 
not be scheduled at the optimal phases. In practice, we have 
a series of RV data points and then the upcoming measure- 
ments are intended to be acquired to yield the smallest un- 
certainty in the orbital eccentricity (or any other parameters 
that are the points of interest) . Since this set of phases yields 
enormous amount of free parameters on which the phases of 
the upcoming measurements are depend, these phases have 
to derived independently by hand for each case. The optsn 
cod^fl is intended to derive these optimal phaseilj] for arbi- 
trary set of fixed observations. 

Since in a multiple planetary system, secular perturba- 
tions result the highest variations in the orbital elements 
(k, h) , the proper selection of RV measurement instants can 
significantly increase the detection probability of other com- 
panions. 



5 SUMMARY 

Although the Newtonian gravitational two-body problem 
(a.k.a., Kepler's problem) is integrable, its solution requires 
transcendent equations. Moreover the parameter of these 
equations has a finite discontinuity in the limit of circu- 
lar orbits, and this ill-behaved property leads to numerical 
disadvantages. For instance the calculation of the optimal 
phases, as presented in Section|31 cannot be performed if the 
orbit was parametrized by (e, -ou) instead of (fc, h) since the 
partial derivative dv/de\e^o, which is required to be known 
to obtain the Fisher- matrix, does not exist. 

In this paper a new framework has been presented to 
describe the time evolution of Kepler's problem using ana- 
lytic functions which omits singular parametrization. This 
formalism is then used to construct orbital solution for three 
dimensional orbits in a similar manner. Using this analytic 
solution, it is straightforward to derive the radial velocity 
function, which is one of the most relevant observable quan- 
tities in the physics of stellar binaries and extrasolar planets. 
The functions p(., ., .) and q(., ., .) can be handled as an ana- 
lytic function, just as simply as if they were a trigonometric 
function in practice. It is also straightforward to carry out 
a Fisher analysis of RV data to estimate the expected un- 
certainty of the physical parameters. The implementation of 
the functions p(., ., .) and q(., ., .) are simple in any kind of 
programming languages, a demonstration code is provided 
in the gnuplot languagqj. 
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* http:/ /szofi. elte.hu/~apal/utils/astro/eof/eof.gnuplot 
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APPENDIX A: THE PARTIAL DERIVATIVES OF THE ECCENTRIC OFFSET FUNCTIONS 



We know that the eccentric ofTsets q and p are continuous functions of A, k and h. Now the partial derivatives of g = q(A, k, h) 
and p = p(A, k, h) are calculated. If all of the partial derivatives are non-singular around zero (where k = h = Qi), one could 
conclude that the eccentric offsets are not only continuous but smooth functions. For the derivation of the partial derivatives 

^ dq dq dq \ 



d{q,p) 
d{X,k,h) 



dp 
V dX 



J 



d\ dk dh 

dp dp 
'dk 'dh 

we will use implicit function theorem. Let us define 

+g cos p + p sin p~ (k cos \ + h sin A) 
— qsinp + pcosp — (k sin X — h cos A) 



Fx(q,p;X,k,h) 

Fy{q,p;\,k,h) 

For a fixed value of A, k and h, Kepler's equation is equivalent with 



Fx{q,p;X,k,h) 
Fy{q,p;X,k,h) 



= 0. 



According to the implicit function theorem, the required partial derivatives are 

dF, dFx ^ 



/ 


dq 


dq 


dq 




dX 


'dk 


'dh 


\ 


dp 


dp 


dp 


dX 


'dk 


'dh 



/ OF^ OF^ \ 
dq dp 



J 



dFy dFy 



V 



dX 



J 



dF. 
dk 



dFy dFy 



dX 



dq dp 

The partial derivatives on the right-hand side of equation 
dFx_ \ 

+ cosp — gsinp + sinp + p cosp 
— sinp + cosp — psinp — gcosp 



dk 
51) are 



dF. 
dh 

dFy_ 
dh 



dq 



dp 



dFy dFy 



dq 
and 

/ ^ 

dX 



dp / 



dF, 
dk 



\ 



dFy dFy 



dFx 
dh 

dFy 



) 



-\-k sin A - 
— fc cos A 



h cos A 
h sin A 



cos A 
- sin A 



— sin A 
+ cos A 



dX dk dh 

Therefore, expanding the above partial derivatives, one can obtain that 

/ dq dq dq \ 
' dX dk dh ^ 



dp 
V 9A 



dp 
'dk 



dp 
'dh / 



1 



+q 



cos(A + p) — k 
+ sin(A -I- p) 



sin(A + p) — h 
— cos(A -l-p) 



(Al) 



(A2) 
(A3) 



(A4) 



(A5) 



(A6) 



(A7) 



(A8) 
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